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Processing seismic data representing a physical system 

The present invention relates to processing seismic data representing a physical system. 
Such processing may be used in data inversion and a particular example of this is the 
5 inversion of time-lapse seismic data. 

Seismic reflection is a technique used to determine details of structures beneath the 
surface of the Earth. The resolution that may be achieved makes this technique the 
method of choice for oil exploration and mapping of subsurface rock structures. It is 
10 also applicable to experimental research that probes the fine structure within the Earth's 
crust and at the crust-mantle boundary. 

The technique involves generating downward-propagating seismic waves in succession 
at a number of locations within the region being explored. A large number of receivers 
1 5 are positioned at intervals away from each source location and these receivers record the 
amplitudes (for example, in terms of pressure, displacement or its derivative) of seismic 
waves reflected back up to the surface from subsurface inhomogeneities over a period 
of time. The recorded waves are usually deconvolved, removing the effects of the 
source and receiver (which have their own response functions). 

20 

Reflection data typically have low amplitudes and are contaminated by multiple 
reflections and other kinds of noise. Various acquisition and processing techniques may 
be used to improve signal-to-noise ratios, such as averaging (stacking) of traces with the 
same midpoint, taking into account different distances between source and receiver, and 

25 discrimination of multiple reflections based on either their periodicity or wavefront 
angles which differ from the primary reflections. Further, the data may be correctly 
positioned in space by a process called migration, which moves dipping events into 
their correct position. When comparisons are made between two or more datasets over 
the same area, careful analysis between the amplitude, time and other attributes of the 

3 0 datasets may be made. 

After the appropriate corrections, which may further include correction for other known 
environmental variables, the data are combined to provide a graphical representation of 
the subsurface inhomogeneities. 
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Seismic reflection data obtained by field experiments are then processed to obtain a 
three dimensional image of subsurface structures as described above. The three 
dimensions refer to the spatial dimensions "illuminated" by the seismic data. The 
5 vertical axis may represent depth or two-way vertical seismic wave travel time. 

t 

The amplitudes of reflected seismic waves are indicative of the subsurface reflection 
strengths, contaminated by noise. The reflection strength depends upon the reflection 
coefficient, which may be defined as a function of the relative contrasts of the elastic 
10 material properties of the subsurface layers. 

The elastic properties of an isotropic, elastic medium are completely described by three 
parameters, for example the two elastic Lam6 parameters and the density. Other 
parameterisations are possible, for example acoustic impedance, shear impedance and 
15 density. A third example is P-wave velocity, S-wave velocity, and density. The 
transformation between different sets of elastic parameters is well defined and 
straightforward. 

In general, the elastic properties vary spatially. In order to explain the relationship 
20 between the elastic properties and the seismic data it may be convenient to imagine the 
subsurface as a stack of geological layers. The layer properties are described by the 
elastic properties of the rocks within the layers while the seismic data are related to the 
contrasts of the layer properties between successive layers. The seismic data are 
therefore suitable for interpreting subsurface layer structures since Ihey image the 
25 boundaries between the layers. 

Seismic inversion is defined herein as the process of transforming (inverting) seismic 
reflection data to elastic material properties, i.e. taking amplitudes (measurements of 
contrasts) and using them to infer physical layer properties. Numerous different seismic 
30 inversion techniques are known. 

Over a period of time, certain types of rock, known as source rocks, will produce 
hydrocarbons. The produced hydrocarbons are then transferred to and stored in rocks 
known as reservoir rocks through various geological processes. During production of 
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hydrocarbons in a subsurface region, the effective elastic material properties of the 
reservoir rocks change with production time, where production time is the fourth 
dimension in seismic 4D analysis. The changes of the effective elastic properties of the 
reservoir rocks may be caused by changes of the pore fluid saturations in the reservoir 
rocks, but also by pressure and temperature changes. Explained by a simple layer-based 
earth model concept, the properties of the reservoir layer are changed during 
production, implying changes in the reflectivity for the upper and lower reservoir 
interfaces. The measurements taken at a further seismic survey are related to the new 
contrasts at the boundaries between adjacent layers. 



Reservoir changes are often inferred from a comparison of the seismic data (e.g. 
amplitudes of seismic waves reflected at interfaces bounding or within the reservoir) for 
different seismic surveys acquired at different stages of the production. A more direct 
interpretation can be based on difference data. Difference data are established by 
15 subtracting two time-separated seismic surveys covering a common part of the earth. 
The difference data, after the proper time-alignment during pre-processing, represent a 
spatial image of the changes of the relative contrasts between the two different 
acquisition times. 

20 For a three dimensional seismic dataset, the classic inversion problem is to estimate the 
elastic material parameters from the three dimensional seismic data. A natural 
extension of 3D inversion to inversion of time-lapse seismic data (4D) is to invert the 
different 3D datasets separately by a known method, and then subtract the results to 
obtain the changes. 

25 

However, the reliability of 4D interpretations is difficult to assess, and are made by 
qualitative assessment. A full consideration of the uncertainties involved is important 
for making an accurate inference of the changes in the reservoir properties between the 
two seismic surveys. The results of such seismic analysis may be important in reservoir 
30 management in that the inferred reservoir properties are used to evaluate, for example, 
new drilling targets and future drainage strategies. 

Seismic inversion provides quantitative estimates of the elastic reservoir properties. 
However, inversion of noisy seismic data is known to be a difficult and ill-posed 
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procedure. An appropriate assessment of the uncertainties in 4D inversion data has not 
previously been possible. 

Commercial time-lapse inversion techniques have become available, but only with brief 
5 descriptions of the methods. Some results have been published (Mesdag et al, 2003, 
Integrated AVO reservoir characterisation and time-lapse analysis of the Widuri field, 
65 th Mtg., Eur., Assn. Expl. Geophys., Extended Abstracts). Such methods apply 
separate inversions of the data with some constraint between the results, e.g. a common 
background model. The time-lapse change is then calculated from the change in 
10 inverted parameters. Sarkar et al, 2003, On the Inversion of time-lapse seismic data, 
73 rd Ann. Internal Mtg.: Soc. Of Expl. Geophys., 1489-1492, mentions inversion of 
seismic differences, but provides no detail of the implementation. None of these 
inversion techniques provide uncertainty bounds on the results. 

1 5 According to a first aspect of the invention, there is provided a method as defined in the 

* ♦ 

appended claim 1 . 

Further aspects and embodiments of the invention are defined in the other appended 
claims. 

20 

It is thus possible to provide a technique which permits improved inversion of seismic 
data representing a physical system. Such a technique may be used to handle errors 
intrinsic to such data and can provide, for example, probability distributions or 
uncertainty bounds on the results of inversion. 

For a better understanding of the present invention and in order to show how the same 
may be carried into effect, preferred embodiments of the invention will now be 
described, by way of example, with reference to the accompanying drawings in which: 

30 Figure 1 is a flow diagram illustrating a method constituting an embodiment of the 
present invention; 

Figures 2a to 2c illustrate a reservoir model; 

Figures 3 a to 3d illustrate modelled seismic data for first and second surveys and their 
associated noise characteristics; 
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Figures 4a to 4c and 5a to 5c illustrate the results of an example of the method of Figure 

i; 

Figure 6 illustrates the uncertainties in the results illustrated in Figures 4a to 4c and 5a 
to 5c; 

5 Figure 7 is a probability map derived from the results of Figures 5a to 5c; and 

Figure 8 is a block schematic diagram of an apparatus for performing the method of 
Figure 1. 

The embodiments of the present invention relate to a method of data inversion that 
10 operates directly on seismic difference data, and in particular to the difference between 
two sets of measured data representing a system in first and second states. In the 
embodiments described herein, the inversion method estimates the changes of the 
elastic material properties of a region of the Earth containing a hydrocarbon reservoir 
due to production or removal of hydrocarbons. The techniques can be based on both 
15 isotropic and anisotropic models, provided that the reflections may be expressed 
linearly. The inversion method gives estimates of the changes in the parameters of the 
model, in addition to the corresponding uncertainty. The solution is obtained by 
combining the information provided by difference data with knowledge obtained prior 
to inversion. The solution is therefore more robust and less vulnerable to instabilities. 
20 Working on difference data rather than inverting the measured data prior to taking the 
difference is advantageous with respect to uncertainty estimation in that it allows a 
correct quantitative statistical treatment of the uncertainties, eliminating the need for 
qualitative interpretation. 

25 To assess the uncertainty of inversion results, the inversion process is cast in a statistical 
setting. The solution of an inverse problem is not limited to a single best-fitting set of 
model parameters, but also characterises the uncertainty of the inversion results. A 
Bayesian setting is chosen for the inversion, although other methods, e.g. least squares, 
can also be used to solve the inversion problem. Li a Bayesian setting it is possible to 

30 combine available prior knowledge with the information contained in the measured 
data. The solution of a Bayesian inverse problem is represented by the posterior 
distribution, which addresses all questions of nonuniqueness and uncertainty. In 
particular, a Gaussian posterior distribution is completely characterised by a posterior 
expectation and a posterior covariance. 
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A model of the earth for the region under consideration is defined, characterised by 
parameters describing the elastic material properties. The embodiments of the present 
invention implement a forward modelling operator, the details of which shall now be 
5 described. 

As described above, an isotropic elastic medium may be fully described using a set of 
three parameters. The embodiments described below adopt the P-wave velocity a(x, t), 
the S-wave velocity p(x, i) and the density p(x, 0 as the parameterisation variables, 
10 where x is the lateral location and t is the two-way vertical seismic wave travel time. 
An alternative parameterisation which may be used takes the acoustic impedance 
Z P =ap, the shear impedance Z s = fip , and the density p as the parameterisation 

variables instead. 

1 5 The weak contrast reflectivity function r for PP reflections can be written as 

where 6 is the reflection angle, a a =(1 + tan 2 6)12, a p ^-A{filaf sin 2 0, and 
a = (! _ 4(/ 0 / a) 2 sin 2 0)12. The values of a a , a p and a p are defined in accordance 
with a known prior background model of the region under consideration. For zero- 
20 incidence reflections, this reduces to 

r(x,/,0) = In Z,(x,r). 
2 ot 

The techniques described herein are with reference to PP reflections, but the inversion is 
equally applicable to other types of reflection, including PS and SS reflections, and in 
25 general to all linear expressions of reflectivity. 

A model parameter vector m can be defined as 

m(x, /) = [In a(x, t), In fi(x, t), In p(x, t)f , 
where T denotes transpose such that m(x,0 is a column vector. A further row vector a 
30 is defined as 

a(x, t } 0) = K (x, r, 0), a p (x, t, 0), a p (x, /, 0)] 
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such that the reflectivity function given above may be defined by the dot product 

d 

r(x,f, 0) = a(x,t, 0) • — m(x,*) . 

ot 

A discrete representation of the reflectivity r in a given time window for a given set of 
5 reflection angles can then be written as a vector 

d 

r(x) = A(x)- — m(x), 

ot 

where A(x) is a sparse matrix defined by a(x,*,0) and m(x) is a discrete representation of 
the model parameters in location x. 

1 0 The seismic data s can be represented by the convolutional model 

s(x, t, 0) = \w(t, 0)r(x, t - r, 0)dr + e(x, t, 0) , 
where w is the seismic wavelet and e is an error term. The wavelet may be angle- 
dependent, but independent of the lateral position x. The wavelet is assumed to be 
stationary within a limited target window. A seismic angle gather at location x can 

1 5 therefore be written, using all of the above, as 

s(x) = WA(x)Dm(x) + e(x), 

where W is a matrix representation of the angle-dependent wavelets, A(x) is the sparse 
matrix defined by a above, and D is merely a differential operator representing partial 
differentiation (with respect to time) of the material parameters. This may be written in 
20 more compact notation by defining the forward modelling operator G(x) = WA(x)D to 
give a description of a seismic angle gather at location x of 

s = Gm + e. 

Seismic data is obtained from the region under consideration on two separate occasions, 
25 shown at step 1 in the flow diagram of Figure 1. Difference data is formed by 
subtracting a first (baseline) seismic dataset s, from a subsequent (repeat) seismic 
dataset s 2 , shown at step 2 of Figure 1. The method described herein may be applied to 
both the difference of pre-stack seismic data, and the difference of fully stacked or 
partially stacked data. The seismic difference data is discrete and is represented on a 
30 seismic grid. A seismic difference vector d is defined as a collection of seismic 
difference data from a set of grid locations, and is represented as d = s 2 - s,. 
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The spatially distributed model parameters are represented on a grid covering the region 
under consideration, using a model parameter vector m containing a collection of the 
model parameters from a set of grid locations, such as those defined above in terms of 
the P and S-wave velocities and density. The model grid should cover the region of 
seismic data collection (the seismic grid), but the two do not need to entirely coincide. 
In most cases, however, coinciding grids would be a preferable choice. The 
embodiments described herein refer to coinciding grids. The step of defining a 
parameterised model is shown at step 3 of Figure 1. 

A model parameter change vector is defined as the difference between the model 
parameter vector m 2 corresponding to the repeat dataset and the model parameter vector 
m, corresponding to the first dataset, that is 5 = m 2 - m,, corresponding to step 4 of 
Figure 1. 

The model parameter change Vector 8 may be defined (for constant x), using the 
parameterisation above, as 



S(x,0 = 



«,(x,r) A(x,r) A (*>')_ 



where the indexes 1 and 2 refer to the first and second seismic data gathers respectively. 
20 For zero-incidence reflections, this reduces to 



8(x,0 = bi 



Z P2 (x,0 



where Z P , , = afii is the acoustic impedance at time i. 

A linear relationship between the model parameter change vector and the seismic data 
25 difference vector can be written in matrix-vector notation as 

d = G6 + e, 

' where G is the forward modelling operator G(x) defined above, and e is the error term 
being related directly to the difference data d. The formulation of this inversion 
problem represents a new approach to the inversion of time-lapse seismic data. 

30 
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The error term e is predominantly a consequence of seismic noise which may be 
characterized with an error expectation vector and an error covariance matrix £ c . The 
error expectation vector is here set to be a zero vector. The error covariance matrix S e 
specifies the variance for each of the elements in the error term vector e, and the 
5 correlations between the different elements therein. The error covariance matrix can be 
estimated from the seismic difference data in regions not influenced by the production 
under consideration, which is not possible unless data are differenced. The error 
expectation and error covariance are determined prior to the inversion. 

10 The knowledge.and uncertainty about the model parameter change vector 8 prior to the 
inversion are characterized via a model parameter change expectation vector p* and a 
model parameter change covariance matrix S 5 . The model parameter change 
expectation vector and the model parameter change covariance matrix are specified 
prior to inversion, and can be determined by analysis of the effects of fluid substitution 

15 and pressure changes due to production. In regions not affected by production, the 
model parameter change expectation is zero (no change is expected). The model change 
parameter covariance matrix specifies the variance for each of the elements ]n(a 2 /a,), 
ln(p 2 /pi), m(p 2 /pi) in the model parameter change vector, and the correlations 
therebetween. 

20 

Explicit analytical expressions are calculated for a solution in which the prior model is 
combined with the information provided by the seismic difference data by forming a 
posterior distribution, shown at step 5 of Figure 1. The solution is represented via an 
updated model parameter change expectation vector, step 6 of Figure 1, and an updated 

25 model parameter change covariance matrix, step 7 of Figure 1. The updated model 
parameter change expectation vector provides an estimator of the change in material 
properties, and thus allows quantitative inferences on the nature of the change to be 
made. The uncertainty is represented via this updated covariance matrix. The analytical 
explicit expressions for the solution provide a computationally fast inversion method 

30 which allows for the assessment of the likelihood of changes in the physical system, 
shown at step 8 of Figure 1 . 
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The statistical properties of the model parameters are assumed to be Gaussian, although 
any other statistical distribution, e.g. the Cauchy distribution, may be used. If each 
single parameter is assumed to be Gaussian defined by an expectation value and a 
variance, then the model parameter change vector 8 is Gaussian, defined by the model 
parameter change expectation vector and the model parameter change covariance 
matrix: 

where n m is the dimension of the model parameter change vector, and and E 5 are as 
defined above. 

The statistical properties of the noise are also assumed to be Gaussian. If each single 
seismic noise sample is assumed to be Gaussian defined by an expectation value and a 
variance, then a vector of noise samples is Gaussian defined by an error expectation 
vector and an error covariance matrix as described above: 

d|5~JV„,(G8,2:,), 

where n d is the dimension of the difference data vector d, and G, 8 and E e are as 
described above. 

These two distribution models therefore imply that the distribution for the seismic 
difference data is also Gaussian: 

d~7/„ rf (Hd,S d ), 

with a difference expectation vector = and difference covariance matrix 
= GE 5 G r + £ e . 

Suppressing the location parameter x, the joint distribution for the model parameter 
change vector 8 and the seismic difference data vector d is then: 



8 
d 



fr 



v 



s 

GE 



s 



Under the described Gaussian assumptions for the prior model and the seismic noise, 
30 the solution is represented via the multi-Gaussian distribution: 



WO 2005/066660 PCT/GB2004/050044 

11 

defined by the new combined expectation and covariance. In this case, the expectation 
vector represents an optimal solution for the 4D inversion problem. In Bayesian 
terminology, this is the posterior distribution with a posterior expectation of 

5 and a posterior covariance of 

4 

The posterior expectation is an optimal estimator for 5 calculated from the difference 

■ 

data d, such that 
10 8 = |i^, 

while the uncertainty is represented by the posterior covariance matrix. 

The posterior estimator Sis an estimator for ln(a 2 /ai), ln(p 2 /Pi) and ln(p 2 /pi). In the 
poststack case, given zero incidence reflections, 5 is an estimator for ln(Z P( 2/Zp,i). 
15 Therefore 8 gives a quantitative estimate of the change in material properties, and the 
corresponding posterior covariance matrix gives its uncertainty. These two quantities 
may then be used to infer changes in production in the reservoir, and other related 
reservoir properties. 

20 The method of the above embodiment will now be illustrated by way of example. 

A reservoir model is defined by acoustic impedance models for the baseline and repeat 
seismic surveys, illustrated in Figures 2a and 2b. The reservoir is originally saturated 
by oil, but during production regions of the reservoir are water-flushed. The drainage is 
25 modelled by an increase of the acoustic impedance. The percentage change in acoustic 
impedance is illustrated in Figure 2c, where the brighter-coloured region indicates the 
region where the physical properties have changed. 

For each of the baseline and repeat models, a seismic forward modelling is performed 
30 by convolution. The wavelet is a Ricker wavelet with a 25Hz centre frequency and 
normalised amplitude. The modelled results S| and s 2 are shown in Figures 3a and 3b 
respectively. A combination of coloured noise ei and white noise e 2 of the form 
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e = s * C) +e 2 was added to the difference data d. The coloured noise represents 
coherent noise, for example a non-repeatability of the source energy between the 
baseline and repeat surveys. The white noise is Gaussian with a variance of a\ . The 
seismic difference was simulated with signal-to-noise ratios of S/N=10 5 and 2, shown in 
5 Figures 3c and 3d respectively. 

To test the inversion method, a stationary prior distribution is assumed, i.e. 

Z P2 (x,r) 



20 



25 



= 0, 



which corresponds to a constant acoustic impedance between surveys. A stationary 
10 covariance function with a\ = 0.039 is also assumed, corresponding to a 0.95 prior 
model interval of ±8% acoustic impedance change, which widely includes the modelled 
increase of about 2-4%. The noise covariance and the wavelet used in the inversion are 
consistent with the model. 

15 The result of inversion of the difference data with a S/N ratio of 10 s is shown in Figures 
4a to 4c, where Figure 4a is the true model, Figure 4b is a waveform representation of 
the seismic difference data, and Figure 4c shows the posterior mean solution. With this 
low noise level the acoustic impedance change is retrieved almost exactly, with low 
uncertainty. 



The result of inversion of the difference data with a S/N ratio of 2 is shown in Figures 
5a to 5c, where Figure 5a is the true model, Figure 5b is a waveform representation of 
the seismic difference data, and Figure 5c shows the posterior mean solution. The 
solution is smoother and has higher uncertainty that that of Figure 4c. 

r 

The uncertainties of the inversion results at CDP 1 1 5 are shown in Figure 6a for the low 
noise case (S/N=10 5 ), and in Figure 6b for the high noise case (S/N=2). The confidence 
region in Figure 6b is much wider than in Figure 6a, as would be expected. 

30 Based on the uncertainty bounds of the inversion results, it is possible to provide 
quantitative probabilities of different reservoir states, e.g. drained, undrained, amount of 
vertical change of a hydrocarbon contact. The probabilities can be used in risk value 
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estimation and commercial decisions. In this example, Figure 7 illustrates a probability 
density map of drained regions (with S/N=2), where darker colouring indicates a higher 
probability of having been drained. This corresponds very well to the modelled case of 
Figure 2c. 



5 



The data inversion methods described above may be embodied in a program for 
controlling a computer to perform the inversion. The program may be stored on a 
storage medium, for example hard or floppy discs, CD or DVD-recordable media or 
flash memory storage products. The program may also be transmitted across a 
10 computer network, for example the Internet or a group of computers connected together 
in a LAN. 

The schematic diagram of Figure 8 illustrates a central processing unit (CPU) 13 
connected to a read-only memory (ROM) 10 and a random access memory (RAM) 12. 

15 The CPU is provided with measured data 14 and model parameters 16 via an 
input/output mechanism 15. The CPU then performs the inversion on the provided data 
in accordance with the instructions provided by the program storage (1 1) (which may be 
• a part of the ROM 10) and provides the output, i.e. the updated model parameters and 
uncertainties 17, via the input/output mechanism 15. The program itself, or any of the 

20 inputs and/or outputs to the system may be provided or transmitted to/from a 
communications network 1 8, which may be, for example, the Internet. 



